home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 16 / CU Amiga Magazine's Super CD-ROM 16 (1997-10-16)(EMAP Images)(GB)[!][issue 1997-11].iso / CUCD / Graphics / Ghostscript / source / gsfemu.c < prev    next >
C/C++ Source or Header  |  1996-12-07  |  18KB  |  723 lines

  1. /* Copyright (C) 1989, 1996 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of Aladdin Ghostscript.
  4.   
  5.   Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  6.   or distributor accepts any responsibility for the consequences of using it,
  7.   or for whether it serves any particular purpose or works at all, unless he
  8.   or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  9.   License (the "License") for full details.
  10.   
  11.   Every copy of Aladdin Ghostscript must include a copy of the License,
  12.   normally in a plain ASCII text file named PUBLIC.  The License grants you
  13.   the right to copy, modify and redistribute Aladdin Ghostscript, but only
  14.   under certain conditions described in the License.  Among other things, the
  15.   License requires that the copyright notice and this notice be preserved on
  16.   all copies.
  17. */
  18.  
  19. /* gsfemu.c */
  20. /* Floating point emulator for gcc */
  21.  
  22. /* We actually only need arch.h + uint and ulong, but because signal.h */
  23. /* may include <sys/types.h>, we have to include std.h to handle (avoid) */
  24. /* redefinition of type names. */
  25. #include "std.h"
  26. #include <signal.h>
  27.  
  28. /*#define TEST*/
  29.  
  30. /*
  31.  * This package is not a fully general IEEE floating point implementation.
  32.  * It implements only one rounding mode (round to nearest) and does not
  33.  * generate or properly handle NaNs.
  34.  *
  35.  * The names and interfaces of the routines in this package were designed to
  36.  * work with gcc.  This explains some peculiarities such as the presence of
  37.  * single-precision arithmetic (forbidden by the ANSI standard) and the
  38.  * omission of unsigned-to-float conversions (which gcc implements with
  39.  * truly grotesque inline code).
  40.  *
  41.  * The following routines do not yet handle denormalized numbers:
  42.  *    addd3 (operands or result)
  43.  *    __muldf3 (operands or result)
  44.  *    __divdf3 (operands or result)
  45.  *    __truncdfsf2 (result)
  46.  *    __extendsfdf2 (operand)
  47.  */
  48.  
  49. /*
  50.  * IEEE single-precision floats have the format:
  51.  *    sign(1) exponent(8) mantissa(23)
  52.  * The exponent is biased by +127.
  53.  * The mantissa is a fraction with an implicit integer part of 1,
  54.  * unless the exponent is zero.
  55.  */
  56. #define fx(ia) (((ia) >> 23) & 0xff)
  57. #define fx_bias 127
  58. /*
  59.  * IEEE double-precision floats have the format:
  60.  *    sign(1) exponent(11) mantissa(52)
  61.  * The exponent is biased by +1023.
  62.  */
  63. #define dx(ld) ((ld[msw] >> 20) & 0x7ff)
  64. #define dx_bias 1023
  65.  
  66. #if arch_is_big_endian
  67. #  define msw 0
  68. #  define lsw 1
  69. #else
  70. #  define msw 1
  71. #  define lsw 0
  72. #endif
  73. /* Double arguments/results */
  74. #define la ((const long *)&a)
  75. #define ula ((const ulong *)&a)
  76. #define lb ((const long *)&b)
  77. #define ulb ((const ulong *)&b)
  78. #define dc (*(const double *)lc)
  79. /* Float arguments/results */
  80. #define ia (*(const long *)&a)
  81. #define ua (*(const ulong *)&a)
  82. #define ib (*(const long *)&b)
  83. #define ub (*(const ulong *)&b)
  84. #define fc (*(const float *)&lc)
  85.  
  86. /* Round a double-length fraction by adding 1 in the lowest bit and */
  87. /* then shifting right by 1. */
  88. #define roundr1(ms, uls)\
  89.   if ( uls == 0xffffffff ) ms++, uls = 0;\
  90.   else uls++;\
  91.   uls = (uls >> 1) + (ms << 31);\
  92.   ms >>= 1
  93.  
  94. /* Extend a non-zero float to a double. */
  95. #define extend(lc, ia)\
  96.   ((lc)[msw] = ((ia) & 0x80000000) + (((ia) & 0x7fffffff) >> 3) + 0x38000000,\
  97.    (lc)[lsw] = (ia) << 29)
  98.  
  99. /* ---------------- Arithmetic ---------------- */
  100.  
  101. /* -------- Add/subtract/negate -------- */
  102.  
  103. double
  104. __negdf2(double a)
  105. {    long lc[2];
  106.     lc[msw] = la[msw] ^ 0x80000000;
  107.     lc[lsw] = la[lsw];
  108.     return dc;
  109. }
  110. float
  111. __negsf2(float a)
  112. {    long lc = ia ^ 0x80000000;
  113.     return fc;
  114. }
  115.  
  116. double
  117. __adddf3(double a, double b)
  118. {    long lc[2];
  119.     int expt = dx(la);
  120.     int shift = expt - dx(lb);
  121.     long sign;
  122.     ulong msa, lsa;
  123.     ulong msb, lsb;
  124.  
  125.     if ( shift < 0 )
  126.       { /* Swap a and b so that expt(a) >= expt(b). */
  127.         double temp = a;
  128.         a = b, b = temp;
  129.         expt += (shift = -shift);
  130.       }
  131.     if ( shift >= 54 )    /* also picks up most cases where b == 0 */
  132.       return a;
  133.     if ( !(lb[msw] & 0x7fffffff) )
  134.       return a;
  135.     sign = la[msw] & 0x80000000;
  136.     msa = (la[msw] & 0xfffff) + 0x100000, lsa = la[lsw];
  137.     msb = (lb[msw] & 0xfffff) + 0x100000, lsb = lb[lsw];
  138.     if ( (la[msw] ^ lb[msw]) >= 0 )
  139.       { /* Adding numbers of the same sign. */
  140.         if ( shift >= 32 )
  141.           lsb = msb, msb = 0, shift -= 32;
  142.         if ( shift )
  143.           { --shift;
  144.             lsb = (lsb >> shift) + (msb << (32 - shift));
  145.         msb >>= shift;
  146.         roundr1(msb, lsb);
  147.           }
  148.         if ( lsb > (ulong)0xffffffff - lsa )
  149.           msa++;
  150.         lsa += lsb;
  151.         msa += msb;
  152.         if ( msa > 0x1fffff )
  153.           { roundr1(msa, lsa);
  154.         /* In principle, we have to worry about exponent */
  155.         /* overflow here, but we don't. */
  156.         ++expt;
  157.           }
  158.       }
  159.     else
  160.       { /* Adding numbers of different signs. */
  161.         if ( shift > 53 )
  162.           return a;        /* b can't affect the result, even rounded */
  163.         if ( shift == 0 && (msb > msa || (msb == msa && lsb >= lsa)) )
  164.           { /* This is the only case where the sign of the result */
  165.         /* differs from the sign of the first operand. */
  166.         sign ^= 0x80000000;
  167.         msa = msb - msa;
  168.         if ( lsb < lsa ) msa--;
  169.         lsa = lsb - lsa;
  170.           }
  171.         else
  172.           { if ( shift >= 33 )
  173.           { lsb = ((msb >> (shift - 32)) + 1) >> 1; /* round */
  174.             msb = 0;
  175.           }
  176.             else if ( shift )
  177.           { lsb = (lsb >> (shift - 1)) + (msb << (33 - shift));
  178.             msb >>= shift - 1;
  179.             roundr1(msb, lsb);
  180.           }
  181.         msa -= msb;
  182.             if ( lsb > lsa ) msa--;
  183.         lsa -= lsb;
  184.           }
  185.         /* Now renormalize the result. */
  186.         /* For the moment, we do this the slow way. */
  187.         if ( !(msa | lsa) )
  188.           return 0;
  189.         while ( msa < 0x100000 )
  190.           { msa = (msa << 1) + (lsa >> 31);
  191.             lsa <<= 1;
  192.         expt -= 1;
  193.           }
  194.         if ( expt <= 0 )
  195.           { /* Underflow.  Return 0 rather than a denorm. */
  196.         lc[msw] = sign;
  197.         lc[lsw] = 0;
  198.         return dc;
  199.           }
  200.       }
  201.     lc[msw] = sign + ((ulong)expt << 20) + (msa & 0xfffff);
  202.     lc[lsw] = lsa;
  203.     return dc;
  204. }
  205. float
  206. __addsf3(float a, float b)
  207. {    return (float)((double)a + (double)b);
  208. }
  209.  
  210. double
  211. __subdf3(double a, double b)
  212. {    long nb[2];
  213.     nb[msw] = lb[msw] ^ 0x80000000;
  214.     nb[lsw] = lb[lsw];
  215.     return a + *(const double *)nb;
  216. }
  217. float
  218. __subsf3(float a, float b)
  219. {    return (float)((double)a - (double)b);
  220. }
  221.  
  222. /* -------- Multiplication -------- */
  223.  
  224. double
  225. __muldf3(double a, double b)
  226. {    long lc[2];
  227.     ulong sign;
  228.     uint H, I, h, i;
  229.     ulong p0, p1, p2;
  230.     ulong expt;
  231.  
  232.     if ( !(la[msw] & 0x7fffffff) || !(lb[msw] & 0x7fffffff) )
  233.       return 0;
  234.     /*
  235.      * We now have to multiply two 53-bit fractions to produce a 53-bit
  236.      * result.  Since the idiots who specified the standard C libraries
  237.      * don't allow us to use the 32 x 32 => 64 multiply that every
  238.      * 32-bit CPU provides, we have to chop these 53-bit numbers up into
  239.      * 14-bit chunks so we can use 32 x 32 => 32 multiplies.
  240.      */
  241. #define chop_ms(ulx, h, i)\
  242.   h = ((ulx[msw] >> 7) & 0x1fff) | 0x2000,\
  243.   i = ((ulx[msw] & 0x7f) << 7) | (ulx[lsw] >> 25)
  244. #define chop_ls(ulx, j, k)\
  245.   j = (ulx[lsw] >> 11) & 0x3fff,\
  246.   k = (ulx[lsw] & 0x7ff) << 3
  247.     chop_ms(ula, H, I);
  248.     chop_ms(ulb, h, i);
  249. #undef chop
  250. #define prod(m, n) ((ulong)(m) * (n)) /* force long multiply */
  251.     p0 = prod(H, h);
  252.     p1 = prod(H, i) + prod(I, h);
  253.     /* If these doubles were produced from floats or ints, */
  254.     /* all the other products will be zero.  It's worth a check. */
  255.     if ( (ula[lsw] | ulb[lsw]) & 0x1ffffff )
  256.       { /* No luck. */
  257.         /* We can almost always avoid computing p5 and p6, */
  258.         /* but right now it's not worth the trouble to check. */
  259.         uint J, K, j, k;
  260.  
  261.         chop_ls(ula, J, K);
  262.         chop_ls(ulb, j, k);
  263.         { ulong p6 = prod(K, k);
  264.           ulong p5 = prod(J, k) + prod(K, j) + (p6 >> 14);
  265.           ulong p4 = prod(I, k) + prod(J, j) + prod(K, i) + (p5 >> 14);
  266.           ulong p3 = prod(H, k) + prod(I, j) + prod(J, i) + prod(K, h) +
  267.         (p4 >> 14);
  268.  
  269.           p2 = prod(H, j) + prod(I, i) + prod(J, h) + (p3 >> 14);
  270.         }
  271.       }
  272.     else
  273.       { p2 = prod(I, i);
  274.       }
  275.     /* Now p0 through p2 hold the result. */
  276.     /****** ASSUME expt IS 32 BITS WIDE. ******/
  277.     expt = (la[msw] & 0x7ff00000) + (lb[msw] & 0x7ff00000) -
  278.       (dx_bias << 20);
  279.     /* Now expt is in the range [-1023..3071] << 20. */
  280.     /* Values outside the range [1..2046] are invalid. */
  281.     p1 += p2 >> 14;
  282.     p0 += p1 >> 14;
  283.     /* Now the 56-bit result consists of p0, the low 14 bits of p1, */
  284.     /* and the low 14 bits of p2. */
  285.     /* p0 can be anywhere between 2^26 and almost 2^28, so we might */
  286.     /* need to adjust the exponent by 1. */
  287.     if ( p0 < 0x8000000 )
  288.       { p0 = (p0 << 1) + ((p1 >> 13) & 1);
  289.         p1 = (p1 << 1) + ((p2 >> 13) & 1);
  290.         p2 <<= 1;
  291.       }
  292.     else
  293.       expt += 0x100000;
  294.     /* The range of expt is now [-1023..3072] << 20. */
  295.     /* Round the mantissa to 53 bits. */
  296.     if ( !((p2 += 4) & 0x3ffc) && !(++p1 & 0x3fff) && ++p0 >= 0x10000000 )
  297.       { p0 >>= 1, p1 = 0x2000;
  298.         /* Check for exponent overflow, just in case. */
  299.         if ( (ulong)expt < 0xc0000000 )
  300.           expt += 0x100000;
  301.       }
  302.     sign = (la[msw] ^ lb[msw]) & 0x80000000;
  303.     if ( expt == 0 )
  304.       { /* Underflow.  Return 0 rather than a denorm. */
  305.         lc[msw] = sign;
  306.         lc[lsw] = 0;
  307.       }
  308.     else if ( (ulong)expt >= 0x7ff00000 )
  309.       { /* Overflow or underflow. */
  310.         if ( (ulong)expt <= 0xc0000000 )
  311.           { /* Overflow. */
  312.         raise(SIGFPE);
  313.         lc[msw] = sign + 0x7ff00000;
  314.         lc[lsw] = 0;
  315.           }
  316.         else
  317.           { /* Underflow.  See above. */
  318.         lc[msw] = sign;
  319.         lc[lsw] = 0;
  320.           }
  321.       }
  322.     else
  323.       { lc[msw] = sign + expt + ((p0 >> 7) & 0xfffff);
  324.         lc[lsw] = (p0 << 25) | ((p1 & 0x3fff) << 11) | ((p2 >> 3) & 0x7ff);
  325.       }
  326.     return dc;
  327. #undef prod
  328. }
  329. float
  330. __mulsf3(float a, float b)
  331. {    long da[2], db[2];
  332.  
  333.     if ( !(ia & 0x7fffffff) || !(ib & 0x7fffffff) )
  334.       return 0;
  335.     extend(da, ia);
  336.     extend(db, ib);
  337.     return (float)(*(const double *)da * *(const double *)db);
  338. }
  339.  
  340. /* -------- Division -------- */
  341.  
  342. double
  343. __divdf3(double a, double b)
  344. {    long lc[2];
  345.     /*
  346.      * Multi-precision division is really, really awful.
  347.      * We generate the result more or less brute force,
  348.      * 11 bits at a time.
  349.      */
  350.     ulong sign = (la[msw] ^ lb[msw]) & 0x80000000;
  351.     ulong msa = (la[msw] & 0xfffff) | 0x100000, lsa = la[lsw];
  352.     ulong msb = (lb[msw] & 0xfffff) | 0x100000, lsb = lb[lsw];
  353.     uint qn[5];
  354.     int i;
  355.     ulong msq, lsq;
  356.     int expt = dx(la) - dx(lb) + dx_bias;
  357.  
  358.     if ( !(lb[msw] & 0x7fffffff) )
  359.       { /* Division by zero. */
  360.         raise(SIGFPE);
  361.         lc[lsw] = 0;
  362.         lc[msw] =
  363.           (la[msw] & 0x7fffffff ?
  364.            sign + 0x7ff00000 /* infinity */ :
  365.            0x7ff80000 /* NaN */);
  366.         return dc;
  367.       }
  368.     if ( !(la[msw] & 0x7fffffff) )
  369.       return 0;
  370.     for ( i = 0; i < 5; ++i )
  371.       { uint q;
  372.         ulong msp, lsp;
  373.  
  374.         msa = (msa << 11) + (lsa >> 21);
  375.         lsa <<= 11;
  376.         q = msa / msb;
  377.         /* We know that 2^20 <= msb < 2^21; q could be too high by 1, */
  378.         /* but it can't be too low. */
  379.         /* Set p = d * q. */
  380.         msp = q * msb;
  381.         lsp = q * (lsb & 0x1fffff);
  382.         { ulong midp = q * (lsb >> 21);
  383.           msp += (midp + (lsp >> 21)) >> 11;
  384.           lsp += midp << 21;
  385.         }
  386.         /* Set a = a - p, i.e., the tentative remainder. */
  387.         if ( msp > msa || (lsp > lsa && msp == msa) )
  388.           { /* q was too large. */
  389.         --q;
  390.         if ( lsb > lsp )
  391.           msp--;
  392.         lsp -= lsb;
  393.         msp -= msb;
  394.           }
  395.         if ( lsp > lsa )
  396.           msp--;
  397.         lsa -= lsp;
  398.         msa -= msp;
  399.         qn[i] = q;
  400.       }
  401.     /* Pack everything together again. */
  402.     msq = (qn[0] << 9) + (qn[1] >> 2);
  403.     lsq = (qn[1] << 30) + (qn[2] << 19) + (qn[3] << 8) + (qn[4] >> 3);
  404.     if ( msq < 0x100000 )
  405.       { msq = (msq << 1) + (lsq >> 31);
  406.         lsq <<= 1;
  407.         expt--;
  408.       }
  409.     /* We should round the quotient, but we don't. */
  410.     if ( expt <= 0 )
  411.       { /* Underflow.  Return 0 rather than a denorm. */
  412.         lc[msw] = sign;
  413.         lc[lsw] = 0;
  414.       }
  415.     else if ( expt >= 0x7ff )
  416.       { /* Overflow. */
  417.         raise(SIGFPE);
  418.         lc[msw] = sign + 0x7ff00000;
  419.         lc[lsw] = 0;
  420.       }
  421.     else
  422.       { lc[msw] = sign + (expt << 20) + (msq & 0xfffff);
  423.         lc[lsw] = lsq;
  424.       }
  425.     return dc;
  426. }
  427. float
  428. __divsf3(float a, float b)
  429. {    return (float)((double)a / (double)b);
  430. }
  431.  
  432. /* ---------------- Comparisons ---------------- */
  433.  
  434. static int
  435. compared2(const long *pa, const long *pb)
  436. {
  437. #define upa ((const ulong *)pa)
  438. #define upb ((const ulong *)pb)
  439.     if ( pa[msw] == pb[msw] )
  440.       { int result = (upa[lsw] < upb[lsw] ? -1 :
  441.               upa[lsw] > upb[lsw] ? 1 : 0);
  442.         return (pa[msw] < 0 ? -result : result);
  443.       }
  444.     if ( (pa[msw] & pb[msw]) < 0 )
  445.       return (pa[msw] < pb[msw] ? 1 : -1);
  446.     /* We have to treat -0 and +0 specially. */
  447.     else if ( !((pa[msw] | pb[msw]) & 0x7fffffff) && !(pa[lsw] | pb[lsw]) )
  448.       return 0;
  449.     else
  450.       return (pa[msw] > pb[msw] ? 1 : -1);
  451. #undef upa
  452. #undef upb
  453. }
  454.  
  455. /* 0 means true */
  456. int
  457. __eqdf2(double a, double b)
  458. {    return compared2(la, lb);
  459. }
  460.  
  461. /* !=0 means true */
  462. int
  463. __nedf2(double a, double b)
  464. {    return compared2(la, lb);
  465. }
  466.  
  467. /* >0 means true */
  468. int
  469. __gtdf2(double a, double b)
  470. {    return compared2(la, lb);
  471. }
  472.  
  473. /* >=0 means true */
  474. int
  475. __gedf2(double a, double b)
  476. {    return compared2(la, lb);
  477. }
  478.  
  479. /* <0 means true */
  480. int
  481. __ltdf2(double a, double b)
  482. {    return compared2(la, lb);
  483. }
  484.  
  485. /* <=0 means true */
  486. int
  487. __ledf2(double a, double b)
  488. {    return compared2(la, lb);
  489. }
  490.  
  491. static int
  492. comparef2(long va, long vb)
  493. {    /* We have to treat -0 and +0 specially. */
  494.     if ( va == vb )
  495.       return 0;
  496.     if ( (va & vb) < 0 )
  497.       return (va < vb ? 1 : -1);
  498.     else if ( !((va | vb) & 0x7fffffff) )
  499.       return 0;
  500.     else
  501.       return (va > vb ? 1 : -1);
  502. }
  503.  
  504. /* See the __xxdf2 functions above for the return values. */
  505. int
  506. __eqsf2(float a, float b)
  507. {    return comparef2(ia, ib);
  508. }
  509. int
  510. __nesf2(float a, float b)
  511. {    return comparef2(ia, ib);
  512. }
  513. int
  514. __gtsf2(float a, float b)
  515. {    return comparef2(ia, ib);
  516. }
  517. int
  518. __gesf2(float a, float b)
  519. {    return comparef2(ia, ib);
  520. }
  521. int
  522. __ltsf2(float a, float b)
  523. {    return comparef2(ia, ib);
  524. }
  525. int
  526. __lesf2(float a, float b)
  527. {    return comparef2(ia, ib);
  528. }
  529.  
  530. /* ---------------- Conversion ---------------- */
  531.  
  532. long
  533. __fixdfsi(double a)
  534. {    /* This routine does check for overflow. */
  535.     long i = (la[msw] & 0xfffff) + 0x100000;
  536.     int expt = dx(la) - dx_bias;
  537.  
  538.     if ( expt < 0 )
  539.       return 0;
  540.     if ( expt <= 20 )
  541.       i >>= 20 - expt;
  542.     else if ( expt >= 31 &&
  543.           (expt > 31 || i != 0x100000 || la[msw] >= 0 ||
  544.            ula[lsw] >= 1L << 21)
  545.         )
  546.       { raise(SIGFPE);
  547.         i = (la[msw] < 0 ? 0x80000000 : 0x7fffffff);
  548.       }
  549.     else
  550.       i = (i << (expt - 20)) + (ula[lsw] >> (52 - expt));
  551.     return (la[msw] < 0 ? -i : i);
  552. }
  553.  
  554. long
  555. __fixsfsi(float a)
  556. {    /* This routine does check for overflow. */
  557.     long i = (ia & 0x7fffff) + 0x800000;
  558.     int expt = fx(ia) - fx_bias;
  559.  
  560.     if ( expt < 0 )
  561.       return 0;
  562.     if ( expt <= 23 )
  563.       i >>= 23 - expt;
  564.     else if ( expt >= 31 && (expt > 31 || i != 0x800000 || ia >= 0) )
  565.       { raise(SIGFPE);
  566.         i = (ia < 0 ? 0x80000000 : 0x7fffffff);
  567.       }
  568.     else
  569.       i <<= expt - 23;
  570.     return (ia < 0 ? -i : i);
  571. }
  572.  
  573. double
  574. __floatsidf(long i)
  575. {    long msc;
  576.     ulong v;
  577.     long lc[2];
  578.  
  579.     if ( i > 0 )
  580.       msc = 0x41e00000 - 0x100000, v = i;
  581.     else if ( i < 0 )
  582.       msc = 0xc1e00000 - 0x100000, v = -i;
  583.     else        /* i == 0 */
  584.       return 0;
  585.     while ( v < 0x01000000 )
  586.       v <<= 8, msc -= 0x00800000;
  587.     if ( v < 0x10000000 )
  588.       v <<= 4, msc -= 0x00400000;
  589.     while ( v < 0x80000000 )
  590.       v <<= 1, msc -= 0x00100000;
  591.     lc[msw] = msc + (v >> 11);
  592.     lc[lsw] = v << 21;
  593.     return dc;
  594. }
  595.  
  596. float
  597. __floatsisf(long i)
  598. {    long lc;
  599.  
  600.     if ( i == 0 )
  601.       lc = 0;
  602.     else
  603.       { ulong v;
  604.         if ( i < 0 )
  605.           lc = 0xcf000000, v = -i;
  606.         else
  607.           lc = 0x4f000000, v = i;
  608.         while ( v < 0x01000000 )
  609.           v <<= 8, lc -= 0x04000000;
  610.         while ( v < 0x80000000 )
  611.           v <<= 1, lc -= 0x00800000;
  612.         v = ((v >> 7) + 1) >> 1;
  613.         if ( v > 0xffffff )
  614.           v >>= 1, lc += 0x00800000;
  615.         lc += v & 0x7fffff;
  616.       }
  617.     return fc;
  618. }
  619.  
  620. float
  621. __truncdfsf2(double a)
  622. {    /* This routine does check for overflow, but it converts */
  623.     /* underflows to zero rather than to a denormalized number. */
  624.     long lc;
  625.  
  626.     if ( (la[msw] & 0x7ff00000) < 0x38100000 )
  627.       lc = la[msw] & 0x80000000;
  628.     else if ( (la[msw] & 0x7ff00000) >= 0x47f00000 )
  629.       { raise(SIGFPE);
  630.         lc = (la[msw] & 0x80000000) + 0x7f800000; /* infinity */
  631.       }
  632.     else
  633.       { lc = (la[msw] & 0xc0000000) +
  634.           ((la[msw] & 0x07ffffff) << 3) +
  635.           (ula[lsw] >> 29);
  636.         /* Round the mantissa.  Simply adding 1 works even if the */
  637.         /* mantissa overflows, because it increments the exponent and */
  638.         /* sets the mantissa to zero! */
  639.         if ( ula[lsw] & 0x10000000 )
  640.           ++lc;
  641.       }
  642.     return fc;
  643. }
  644.  
  645. double
  646. __extendsfdf2(float a)
  647. {    long lc[2];
  648.  
  649.     if ( !(ia & 0x7fffffff) )
  650.       lc[msw] = ia, lc[lsw] = 0;
  651.     else
  652.       extend(lc, ia);
  653.     return dc;
  654. }
  655.  
  656. /* ---------------- Test program ---------------- */
  657.  
  658. #ifdef TEST
  659.  
  660. #include <stdio.h>
  661. #include <stdlib.h>
  662.  
  663. int
  664. test(double v1)
  665. {    double v3 = v1 * 3;
  666.     double vh = v1 / 2;
  667.     double vd = v3 - vh;
  668.     double vdn = v1 - v3;
  669.  
  670.     printf("%g=1 %g=3 %g=0.5 %g=2.5 %g=-2\n", v1, v3, vh, vd, vdn);
  671.     return 0;
  672. }
  673.  
  674. float
  675. randf(void)
  676. {    int v = rand();
  677.     v = (v << 16) ^ rand();
  678.     if ( !(v & 0x7f800000) )
  679.       return 0;
  680.     if ( (v & 0x7f800000) == 0x7f800000 )
  681.       return randf();
  682.     return *(float *)&v;
  683. }
  684.  
  685. int
  686. main(int argc, char *argv[])
  687. {    int i;
  688.     test(1.0);
  689.     for ( i = 0; i < 10; ++i )
  690.       { float a = randf(), b = randf(), r;
  691.         int c;
  692.         switch ( (rand() >> 12) & 3 )
  693.           {
  694.           case 0: r = a + b; c = '+'; break;
  695.           case 1: r = a - b; c = '-'; break;
  696.           case 2: r = a * b; c = '*'; break;
  697.           case 3: if ( b == 0 ) continue; r = a / b; c = '/'; break;
  698.           }
  699.         printf("0x%08x %c 0x%08x = 0x%08x\n",
  700.            *(int *)&a, c, *(int *)&b, *(int *)&r);
  701.       }
  702. }
  703.  
  704. /*
  705.    Results from compiling with hardware floating point:
  706.  
  707. 1=1 3=3 0.5=0.5 2.5=2.5 -2=-2
  708. 0x6f409924 - 0x01faa67a = 0x6f409924
  709. 0x00000000 + 0x75418107 = 0x75418107
  710. 0xe32fab71 - 0xc88f7816 = 0xe32fab71
  711. 0x94809067 * 0x84ddaeee = 0x00000000
  712. 0x2b0a5b7d + 0x38f70f50 = 0x38f70f50
  713. 0xa5efcef3 - 0xc5dc1a2c = 0x45dc1a2c
  714. 0xe7124521 * 0x3f4206d2 = 0xe6ddb891
  715. 0x8d175f17 * 0x2ed2c1c0 = 0x80007c9f
  716. 0x419e7a6d / 0xf2f95a35 = 0x8e22b404
  717. 0xe0b2f48f * 0xc72793fc = 0x686a49f8
  718.  
  719.  */
  720.  
  721.  
  722. #endif
  723.